function [IA_Network, consistent, iterations] = IA_Network_Path_Consistency_Matrix2(IA_Network, i, j, IA_Relation)
%IA_Network_Path_Consistency To Add R(i, j)
%   Detailed explanation goes here

% Todo (RF): both 'N' and 'R' matrices are at least 2x bigger than required 
% as the relationship matrix is symmetric around the (top left to bottom 
% right) diagonal. The diagonal also contains no information as each 
% interval can only be equal to itself and no other relation is possible.

consistent = 1;
iterations = 0;

if IA_Relation_Is_Dummy(IA_Relation)
  return;
end

if IA_Relation_Is_Null(IA_Relation)
  consistent = 0;
  return;
end

N = IA_Network.relations;
N(i, j) = IA_Relation_Intersect(N(i, j), IA_Relation);
N(j, i) = IA_Relation_Invert(N(i, j));

if IA_Relation_Is_Null(N(i, j))
  consistent = 0;
  return;
end

done = 0;
while ~done
  % N^2 = N o N
  parNoN(1:size(N,1), 1:size(N,1), 1:size(N,2)) = IA_Relation_Dummy();
  for k = 1:size(N,1)
    for i = 1:size(N,1)
      for j = 1:size(N,2)
        if (IA_Relation_Has_Less(N(i, k)) && IA_Relation_Has_Greater(N(k, j)))...
        || (IA_Relation_Has_Greater(N(i, k)) && IA_Relation_Has_Less(N(k, j)))...
        || (IA_Relation_Has_During(N(i, k)) && IA_Relation_Has_Contains(N(k, j)))
          continue;
        end
        
        parNoN(k, i, j) = IA_Relation_Constraints_Hogge(N(i, k), N(k, j), N(i, j)); %#ok<AGROW>
        if IA_Relation_Is_Null(parNoN(k, i, j))
          consistent = 0;
          return;
        end
      end
    end
  end
  
  comNoN(1:size(N,1), 1:size(N,2)) = IA_Relation_Dummy();
  for k = 1:size(N,1)
    for i = 1:size(N,1)
      for j = 1:size(N,2)
        comNoN(i, j) = IA_Relation_Intersect(comNoN(i, j), parNoN(k, i, j)); %#ok<AGROW>
      end
    end
  end
  NoN = comNoN;
  
  done = 1;
  for i = 1:size(N,1)
    for j = 1:size(N,2)
      % IF N == N ^ 2 THEN QUIT
      if ~IA_Relation_Is_Equal(N(i, j), NoN(i, j))
        done = 0;
      end
      
      % N = N . N ^ 2
      N(i, j) = IA_Relation_Intersect(N(i, j), NoN(i, j));
      if IA_Relation_Is_Null(N(i, j))
        consistent = 0;
        return;
      end
    end
  end
  count = count + 1;
end

IA_Network.relations = N;
end